#!/bin/csh
# overerplot_timeseries.csh
##############################################################
# Edit overplot-timeseries_files
# Run the script using: ./overplot_timeseries.csh time0 time1 field
# where time0 and time1 are initial and final times of the series
# and field is the field to plot (temperature, nusselt, velocity, mixfrac)
# NB: the ouput files are stored in the dir1 directory selected in
#     overplot-timeseries_files under ../Plots/Comparisons/
##########################################################################

#---------------------------
# Options
#---------------------------
set makejpg = 1
set jpgdensity = 200
set deletepdf = 1
set dimx = 20
set dimz = 10

#-------------------------------------------------------------------------
# Read file containing the simulation directories and the profiles to plot
#-------------------------------------------------------------------------
source overplot-timeseries_files

#------------------------------
# Set ouput directories
#------------------------------
set dirplots0 = {$dir1}Plots/
set dirplots = {$dir1}Plots/Comparisons/ 

#------------------------------
# Set timeseries files
#------------------------------
set f1 = ${dir1}timeseries.dat
set f2 = ${dir2}timeseries.dat
set f3 = ${dir3}timeseries.dat
set f4 = ${dir4}timeseries.dat
set f5 = ${dir5}timeseries.dat
set f6 = ${dir6}timeseries.dat
set f7 = ${dir7}timeseries.dat
set f8 = ${dir8}timeseries.dat
set f9 = ${dir9}timeseries.dat
set f10 = ${dir10}timeseries.dat

#---------------------------------------------------------
# Prepare arrays of timeseries files and gmt plot options 
#---------------------------------------------------------
set ts = ($f1 $f2 $f3 $f4 $f5 $f6 $f7 $f8 $f9 $f10)
set opt = ($opt1 $opt2 $opt3 $opt4 $opt5 $opt6 $opt7 $opt8 $opt9 $opt10)

#---------------------------
# Create output directory
#---------------------------
if (! -d $dirplots0) then
  mkdir $dirplots0
endif

if (! -d $dirplots) then
  mkdir $dirplots
endif

#---------------------------
# Select field to plot
#---------------------------
set field = ${3}
if (${3} !~ temperature && ${3} !~ nusselt && ${3} !~ velocity && ${3} !~ mixfrac) then
  echo 'The field you want to plot is either not (yet) supported in this routine or has been misspelled'
  exit
endif

#----------------------------
# Set plot origin (in cm)
#----------------------------
set X0 = 3
set Y0 = 5

#------------------------
# Set min and max time
#------------------------
set time0 = ${1}
set time1 = ${2}
set atime = `echo "scale=3; ($time1 - $time0)/5" | bc -l`

set output = t.ps

#----------------------------------
# Set gmtdefaults
# ---------------------------------
source setgmtdefaults

####################
# Temperature
####################
if ($field == temperature) then
  echo "Overplot of temperature time-series of..." 

  set pdfoutput = temperature_ts_overplot.pdf
  set jpgoutput = temperature_ts_overplot.jpg

  set title = ""
  set xlabel = "Time"
  set zlabel = "Temperature"

  set min = 0
  set max = 1
  set delta = 0.2
endif

####################
# Nusselt
####################
if ($field == nusselt) then
  echo "Overplot of Nusselt number time-series of..." 

  set pdfoutput = nusselt_ts_overplot.pdf
  set jpgoutput = nusselt_ts_overplot.jpg

  set title = ""
  set xlabel = "Time"
  set zlabel = "Nusselt@~\040@~number"

  set min = 0
  set max = 30
  set delta = 5
endif

####################
# Velocity
####################
if ($field == velocity) then
  echo "Overplot of RMS velocity time-series of..." 

  set pdfoutput = velocity_ts_overplot.pdf
  set jpgoutput = velocity_ts_overplot.jpg

  set title = ""
  set xlabel = "Time"
  set zlabel = "Velocity"

  set min = 0
  set max = 1500
  set delta = 250
endif

####################
# Mixed fraction
####################
if ($field == mixfrac) then
  echo "Overplot of mixed fraction time-series of..." 

  set pdfoutput = mixfrac_ts_overplot.pdf
  set jpgoutput = mixfrac_ts_overplot.jpg

  set title = ""
  set xlabel = "Time"
  set zlabel = "Mixed@~\040@~fraction"

  set min = 0
  set max = 1
  set delta = 0.2
endif

#-------------------------
# GMT plot
#-------------------------
psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba${atime}:${xlabel}:/a${delta}:${zlabel}:WSen:.${title}: -K >! $output
  
set i = 1
while ( $i <= $nts )
  echo $ts[$i]
  
  set K = -K
  if ($i == $nts) then
    set K = ""
  endif

  if ($field == temperature) then
    awk 'NR>1 {print $1, $3}' $ts[$i]  >! p
  else if ($field == nusselt) then
    awk 'NR>1 {print $1, $2}' $ts[$i]  >! p
  else if ($field == velocity) then
    awk 'NR>1 {print $1, $4}' $ts[$i]  >! p
  else if ($field == mixfrac) then
    awk 'NR>1 {print $1, $16}' $ts[$i]  >! p
  endif

  psxy p -R -JX $opt[$i] $K -O >> $output

  @ i = $i + 1
end


#----------------------------------------
# Conversions and cleaning
#----------------------------------------
ps2raster $output -Au -P -Tf -F$pdfoutput
if ($makejpg == 1) then
  ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  mv -f *.jpg $dirplots
endif

if ($makejpg == 1 && $deletepdf == 1) then
  rm -f *.pdf
else 
  mv -f *.pdf $dirplots
endif

rm -f p *.ps 
